suppressPackageStartupMessages(library(tidyverse))
library(readstata13)
library(BAS)

source("Code/silent.R")

silent(
    data_electoral_turnover <- haven::read_dta("Data/Covid19_Incumbency_Dataset_workfile_25Nov24.dta") %>%
        as_tibble()
)

silent(
    lm_fit_all_controls_BAS <- bas.lm(
        formula = excessmort_vb_remain ~ turnover + log_MaleToFemaleRatio + log_shareImmigrants + log_shareBlueCollar + log_PopDensity + log_sq_PopDensity + log65 + log75 + log80 + log6575 + log6580 + log7580 + log3 + log_MedianStandardLiving10k + log_basemort_vb_remain + log_shareUnemp + delta_Immigrants + delta_Unemp + delta_logStandardLiving + log_nl2019 + delta_log_nl + ncandidates,
        data = data_electoral_turnover,
        prior = "JZS",
        modelprior = uniform(),
        initprobs = "eplogp",
        force.heredity = FALSE, pivot = TRUE, include.always = ~ 1 + turnover
    )
)
lm_fit_all_controls_BAS %>%
    image(subset = 1:23)

pdf("Figures/figA9.pdf")
lm_fit_all_controls_BAS %>%
    image(rotate = FALSE, subset = 1:23)
dev.off()

coef.lm_fit_all_controls_BAS <- coef(lm_fit_all_controls_BAS)

pdf("Figures/figA10a.pdf")
plot(coef.lm_fit_all_controls_BAS, subset = 2, ask = F)
dev.off()

pdf("Figures/figA10b.pdf")
plot(coef.lm_fit_all_controls_BAS, subset = 3, ask = F)
dev.off()
